Structure Analysis

Examined "global" ancestry for all available samples. Old replicates and some suspected contamination/switching have been removed. The masterpopmap has ALL of these samples.

Admixture model, LD=0, allele freqs are uncorrelated, alpha is inferred for each popualtion.

Read in Results and set up popmap/poporder for visualization and downstream analysis.

sfiles <- list.files(path="Results/",
  pattern=glob2rx("outfile_*f"),full.names=TRUE)

slist <- readQ(sfiles,filetype="structure",
               indlabfromfile = TRUE,
               readci = TRUE)

popmap <- left_join(data.frame(sampID = rownames(slist[[1]])),
                    masterpopmap)
## Joining, by = "sampID"
popID <- data.frame(popID = popmap$popID)
popID <- data.frame(lapply(popID, as.character), stringsAsFactors=FALSE)

pop_order <- sample_meta.dat[,c("popID", "Order")]
pop_order <- arrange(pop_order, Order)[,1]

Sample size for just Structure analysis

samplesize <- left_join(popmap, sample_meta.dat) %>%
  count(popID) %>%
  merge(sample_meta.dat)
## Joining, by = "popID"
samplesize
##         popID  n
## 1        11CO 12
## 2        12AZ 15
## 3        15TX 16
## 4        16TX 18
## 5        18TX 17
## 6        19TX 15
## 7         1WY 17
## 8        20NM 17
## 9        21NM 14
## 10       26NM 18
## 11       27NM 20
## 12       28NM 17
## 13        2CO 19
## 14       31NV  7
## 15       32UT 20
## 16       33NV  9
## 17       34UT 16
## 18       37CR  5
## 19       38CR 11
## 20       39CR 12
## 21       41CR 16
## 22       43GR 15
## 23       44GR 16
## 24       46CH 24
## 25       47TX 25
## 26       49TX 27
## 27        4CO 11
## 28       50TX  8
## 29       51TX 14
## 30       52TX 21
## 31        5CO 19
## 32        6KS 10
## 33        8OK 20
## 34 CARINA_LAB 14
## 35    SUB_LAB 17
##                                                               Name Latitude
## 1                                                        Wilkinson    37.34
## 2                                                         Big Bend    35.12
## 3                                                        NE Post\n    33.32
## 4                                                   Lake JB Thomas    32.61
## 5                                                             Orla    31.49
## 6                                                        Aspermont    33.17
## 7                                                           Lovell    44.86
## 8                                                         Malaga\n    32.22
## 9                                         Artesia Wildlife Reserve    32.98
## 10                                                     Lake Sumner    34.15
## 11                                                       Roswell E    33.40
## 12                                                  Tucumcari Lake    35.19
## 13                                                  Fountain Creek    38.34
## 14                                       Virgin River - Gold Butte    36.69
## 15                                                      St. George    37.07
## 16                                          Lake Mead, Stewarts Pt    36.38
## 17                                                           Delta    39.23
## 18                                         Rethimno, Crete, Greece    35.37
## 19                                                       Crete, GR    35.42
## 20                                       Panaramnos, Crete, Greece    35.42
## 21                                          Plakias, Crete, Greece    35.19
## 22                                                  Posidi, Greece    39.97
## 23                                            Delta Aksiou, Greece    40.55
## 24                                                    Bitun, China    47.30
## 25 Rio Grande Village - from T. cheninsis x T. ramosissima hybrids    29.18
## 26                                                     Santa Elana    29.16
## 27                                   Adobe Reservoir (Blue Lake)\n    38.26
## 28                                                    Presidio Hwy    29.34
## 29                                             CO Canyon Boat Ramp    29.34
## 30                                                        Tornillo    31.44
## 31                                                       SE corner    37.70
## 32                                                 W Finney County    37.99
## 33                                                          Guymon    36.70
## 34                                       D. carinata Lab Culture\n       NA
## 35                                     D. sublineata Lab Culture\n       NA
##    Longitude Date.Collected Order                   Grp              cat
## 1    -104.16        9/18/14     7            carinulata        expansion
## 2    -114.64        9/18/14    11            carinulata        expansion
## 3    -101.26        9/26/14    27 potential_hybrid_zone        expansion
## 4    -101.22        9/26/14    30 potential_hybrid_zone original_release
## 5    -103.48        9/27/14    32 potential_hybrid_zone original_release
## 6    -100.24        9/27/14    28 potential_hybrid_zone        expansion
## 7    -108.18         9/3/14     2            carinulata original_release
## 8    -104.08        9/27/14    31 potential_hybrid_zone original_release
## 9    -104.44        10/7/14    29 potential_hybrid_zone original_release
## 10   -104.48        9/29/14    25 potential_hybrid_zone original_release
## 11   -104.41        9/29/14    26 potential_hybrid_zone original_release
## 12   -103.69        10/6/14    24 potential_hybrid_zone original_release
## 13   -104.61         9/2/14     4            carinulata original_release
## 14   -114.26       10/15/14     9            carinulata        expansion
## 15   -113.58       10/15/14     8            carinulata        expansion
## 16   -114.40       10/15/14    10            carinulata        expansion
## 17   -112.93       10/15/14     3            carinulata original_release
## 18     24.47         7/2/15    18       native_elongata           native
## 19     24.69         7/4/15    17       native_elongata           native
## 20     24.68         7/4/15    16       native_elongata           native
## 21     24.40         7/4/15    20       native_elongata           native
## 22     23.37         7/6/15    15       native_elongata           native
## 23     22.74         7/6/15    14       native_elongata           native
## 24     87.75         7/3/16     1            carinulata           native
## 25   -102.96         8/3/17    36 potential_hybrid_zone original_release
## 26   -103.60         8/3/17    37 potential_hybrid_zone original_release
## 27   -103.25         9/4/14     5            carinulata        expansion
## 28   -104.07         8/3/17    34 potential_hybrid_zone original_release
## 29   -104.06         8/3/17    35 potential_hybrid_zone original_release
## 30   -106.09         8/3/17    33 potential_hybrid_zone        expansion
## 31   -103.42        8/21/14     6            carinulata        expansion
## 32   -101.08         8/8/14    22 potential_hybrid_zone        expansion
## 33   -101.55        9/16/14    23 potential_hybrid_zone        expansion
## 34        NA                   12                   lab              lab
## 35        NA                   13                   lab              lab
write.table(samplesize, "pop_metadata_samplesize.tsv", sep = "\t", row.names = F, quote = F)

What is the best K?

sr1 <- summariseQ(tabulateQ(slist))
p <- evannoMethodStructure(data=sr1,exportplot=F,returnplot=T,returndata=F,basesize=12,linesize=0.7)
grid.arrange(p)

evannoMethodStructure(data=sr1,
                  exportplot=T,
                  exportpath = "plots",
                outputfilename = "global_deltaK",
                returndata=F,
                height = 10, width = 9)
## plots/global_deltaK.png exported.

K=2 and K=3 are super stable. D. carinata appears as a mix between sublineata and elongata. I am surprised to see, for the first time in these samples, indications of admixture between carinlata and sublineata along the Rio Grande.

As usual, K=4 reflects lab populations, elongata in the native range, and hints at carinulata substructure.

k4 <- plotQ(alignK(slist[c(41:50)]),
            imgoutput="join",
            returnplot=T,exportplot=F,basesize=11,
            clustercol=cbbPalette,
            grplab = popID,
            ordergrp = T,
            grplabangle=-90)
grid.arrange(k4$plot[[1]])

Code to save a nice version of all reps of K=4 for the manuscript.

fn1 <- function(x) attr(x,"k")
spnames <- rep("", 10)

plotQ(alignK(slist[c(41:50)]),
            imgoutput="join",
            exportpath = "plots",
  outputfilename =  "Joined-K4",  
  exportplot=T,
  showyaxis=F,
  showsp = T,
  splab=rep("", 10),
  clustercol=cbbPalette,
  subsetgrp = pop_order,
  height = .005,
  ordergrp = T,
  grplab = popID,
  linepos = 1,
  divtype = 1,
  grplabangle=90,
  grplabpos = .8,
  grplabspacer = 0,
  basesize=14,
  grplabsize = 15,
  width = 150,
  grplabheight=150)
## Drawing plot ...
## plots/Joined-K4.png exported.

Assign species identity to each cluster

Extract one q file that reflects known species identity

global_k4_q <- slist[[41]]
global_k4_q$sampID <- row.names(global_k4_q)

## which cluster is carinulata?
gr_samples <- masterpopmap[grep( "CR", masterpopmap[,2]), ]
gr_q <- semi_join(global_k4_q, gr_samples)
max_cl_gr <- sapply(gr_q[,1:4], max)
elongata_cluster <- which(colnames(gr_q) == colnames(gr_q[which(max_cl_gr==1)]))
elongata_cluster
## [1] 1
## which cluster is elongata?
ch_samples <- masterpopmap[grep( "CH", masterpopmap[,2]), ]
ch_q <- semi_join(global_k4_q, ch_samples)
max_cl_ch <- sapply(ch_q[,1:4], max)
carinulata_cluster <- which(colnames(ch_q) == colnames(ch_q[which(max_cl_ch==1)]))
carinulata_cluster
## [1] 3
## which cluster is sublineata?
sub_samples <- masterpopmap[grep( "SUB_LAB", masterpopmap[,2]), ]
sub_q <- semi_join(global_k4_q, sub_samples)
max_cl_sub <- sapply(sub_q[,1:4], max)
sublineata_cluster <- which(colnames(sub_q) == colnames(sub_q[which(max_cl_sub==1)]))

## which cluster is carinata?
carina_samples <- masterpopmap[grep( "CARINA_LAB", masterpopmap[,2]),]
carina_q <- semi_join(global_k4_q, carina_samples)
max_cl_cara <- sapply(carina_q[,1:4], max)
carinata_cluster <- which(colnames(carina_q) == colnames(carina_q[which(max_cl_cara==1)]))

cls <- c(carinata_cluster, carinulata_cluster, sublineata_cluster, elongata_cluster)
spp <- c("carinata", "carinulata", "sublineata", "elongata")
clusters <- as.data.frame(cbind(spp = spp[order(cls)], cluster = seq(1,4)))

colnames(global_k4_q)[1:4] <- as.character(clusters$spp)
q_popmap <- merge(global_k4_q, popmap)
head(q_popmap)
##    sampID elongata sublineata carinulata carinata popID
## 1      D1    1.000      0.000       0.00    0.000  38CR
## 2     D12    1.000      0.000       0.00    0.000  41CR
## 3     D13    1.000      0.000       0.00    0.000  41CR
## 4 D15_rep    0.970      0.000       0.03    0.000  41CR
## 5     D16    1.000      0.000       0.00    0.000  41CR
## 6     D17    0.998      0.001       0.00    0.001  43GR
legend <- merge(legend_colors, clusters)
clustercols <- arrange(legend, cluster) %>%
  select(cbbPalette)
clustercols <- clustercols$cbbPalette

legendlab <- arrange(legend, cluster) %>%
  mutate(label = paste("D. ", spp))
legendlab <- legendlab$label[1:4]

Code to save a nice version of one rep of K=4 for the manuscript.

plotQ(slist[41],
            exportpath = "plots",
        outputfilename =  "MainTextFigure_K4",  
  exportplot=T,
  showyaxis=T,
  showsp = F,
  clustercol=clustercols,
  subsetgrp = pop_order,
  sortind="all",
  grplab = popID,
  linepos = 1,
  linesize=0.8,
  pointsize = 5,
  divtype = 1,
  divsize=1.1,
  grplabangle= 90,
  grplabpos = .9,
  grplabspacer = 0,
  grplabheight=35,
  grplabsize = 10,
  height = .05,
  width=100,
  showlegend = T, legendkeysize=30,legendtextsize=30, legendpos="left", legendlab=legendlab,
  # titlelab="Global Population Structure K=4",
  # titlesize = 20,
  basesize=30)

Code to save a nice version of K=2-4 for the manuscript.

# custom strip panel label showing k only
fn1 <- function(x) attr(x,"k")
spnames <- paste0("K=",sapply(slist,fn1))

# Combine K=2, K=3, K=4
k2k3k4 <- c(21,31,41)
plotQ(slist[k2k3k4],
            imgoutput="join",
            exportpath = "plots",
  outputfilename =  "Joined-K2-K3-K4",  
  exportplot=T,
  showyaxis=F,
  clustercol=clustercols,
  subsetgrp = pop_order,
  ordergrp = T,
  grplab = popID,
  linepos = 1,
  divtype = 1,
  # divsize = 1.2,
  grplabangle= 90,
  grplabpos = .8,
  grplabspacer = 0,
  basesize=14,
  showsp = F,
  grplabsize = 6,
  width = 30,
  grplabheight=40)

Exploratory visualizations that are not presented in the manuscript

Multi-line K=4 to visulize individuals

plotQMultiline(slist[41],
      showticks=T,
      clustercol=cbbPalette,
      grplab = popID,
      ordergrp=TRUE,
      spl = 25,
      showindlab = T,
      showyaxis = T,
      barsize = .9,
      subsetgrp = pop_order,
      basesize=11,
      exportpath = "plots",
      outputfilename = "multilineK4",
      useindlab = T)
## Drawing plot ...
## plots/multilineK4-1.png exported.
## Drawing plot ...
## plots/multilineK4-2.png exported.
## Drawing plot ...
## plots/multilineK4-3.png exported.

All reps of each run to visualize results across K values

k2 <- plotQ(alignK(slist[c(21:30)]),
            imgoutput="join",
            returnplot=T,exportplot=F,basesize=11,
            clustercol=cbbPalette,
            subsetgrp = pop_order,
            ordergrp = T,
            grplab = popID,
            linepos = 1,
            grplabangle=-45,
            grplabpos = 0,
            grplabheight=20)
grid.arrange(k2$plot[[1]])

k3 <- plotQ(alignK(slist[c(31:40)]),
            imgoutput="join",
            returnplot=T,exportplot=F,basesize=11,
            clustercol=cbbPalette,
            ordergrp = T,
            grplab = popID,
            grplabangle=-90)
grid.arrange(k3$plot[[1]])

k5 <- plotQ(alignK(slist[c(51:60)]),
            imgoutput="join",
            returnplot=T,exportplot=F,basesize=11,
            clustercol=cbbPalette,
            grplab = popID,
            ordergrp = T,
            grplabangle=-90)
grid.arrange(k5$plot[[1]])

k6 <- plotQ(alignK(slist[c(61:70)]),
            imgoutput="join",
            returnplot=T,exportplot=F,basesize=11,
            clustercol=cbbPalette,
            grplab = popID,
            ordergrp = T,
            grplabangle=-90)
grid.arrange(k6$plot[[1]])

k7 <- plotQ(alignK(slist[c(71:80)]),
            imgoutput="join",
            returnplot=T,exportplot=F,basesize=11,
            clustercol=cbbPalette,
            grplab = popID,
            grplabangle=-90)
grid.arrange(k7$plot[[1]])

k8 <- plotQ(alignK(slist[c(81:90)]),
            imgoutput="join",
            returnplot=T,exportplot=F,basesize=11,
            # clustercol=cbbPalette,
            grplab = popID,
            grplabangle=-90)
grid.arrange(k8$plot[[1]])

Among individuals, how many pure? How many hybrids, and of what type?

Extract confidence intervals.

ci <- data.frame(attributes(slist[[41]])$ci)
colnames(ci) <- paste0(rep(colnames(global_k4_q)[1:4], each=2),
       "_", rep(c("L", "H"), times = 4))
ci <- cbind(ci, global_k4_q)
ci <- merge(ci, masterpopmap)
head(ci)
##    sampID elongata_L elongata_H sublineata_L sublineata_H carinulata_L
## 1      D1      0.999      1.000            0        0.001        0.000
## 2     D12      0.998      1.000            0        0.001        0.000
## 3     D13      0.998      1.000            0        0.001        0.000
## 4 D15_rep      0.953      0.985            0        0.002        0.015
## 5     D16      0.999      1.000            0        0.001        0.000
## 6     D17      0.989      1.000            0        0.005        0.000
##   carinulata_H carinata_L carinata_H elongata sublineata carinulata carinata
## 1        0.001          0      0.001    1.000      0.000       0.00    0.000
## 2        0.001          0      0.001    1.000      0.000       0.00    0.000
## 3        0.001          0      0.001    1.000      0.000       0.00    0.000
## 4        0.046          0      0.001    0.970      0.000       0.03    0.000
## 5        0.001          0      0.001    1.000      0.000       0.00    0.000
## 6        0.001          0      0.005    0.998      0.001       0.00    0.001
##   popID
## 1  38CR
## 2  41CR
## 3  41CR
## 4  41CR
## 5  41CR
## 6  43GR
write.csv(ci,"ancestry_confidenceintervals.csv",
          quote = FALSE)

Use "diagnostic" samples to define threshold.

diag_sublineata <- c("SUB_LAB", "47TX", "49TX", "50TX", "51TX", "52TX")
diagnostic <- data.frame(popID = diag_sublineata, spp=rep("sublineata", length(diag_sublineata)))
diag_elongata <- c("43GR", "44GR", "37CR", "38CR", "39CR", "41CR")
diagnostic <- rbind(diagnostic, 
      data.frame(popID = diag_elongata, spp=rep("elongata", length(diag_elongata))))

ci %>%
  semi_join(diagnostic) %>%
  arrange(desc(carinulata_L))
## Joining, by = "popID"
##      sampID elongata_L elongata_H sublineata_L sublineata_H carinulata_L
## 1        N7      0.000      0.001        0.888        0.934        0.066
## 2   D24_rep      0.910      0.951        0.000        0.001        0.049
## 3   D15_rep      0.953      0.985        0.000        0.002        0.015
## 4        D1      0.999      1.000        0.000        0.001        0.000
## 5       D12      0.998      1.000        0.000        0.001        0.000
## 6       D13      0.998      1.000        0.000        0.001        0.000
## 7       D16      0.999      1.000        0.000        0.001        0.000
## 8       D17      0.989      1.000        0.000        0.005        0.000
## 9       D18      0.998      1.000        0.000        0.001        0.000
## 10      D19      0.998      1.000        0.000        0.001        0.000
## 11       D2      0.999      1.000        0.000        0.001        0.000
## 12      D20      0.998      1.000        0.000        0.001        0.000
## 13      D21      0.997      1.000        0.000        0.001        0.000
## 14      D22      0.998      1.000        0.000        0.001        0.000
## 15      D23      0.998      1.000        0.000        0.001        0.000
## 16      D27      0.998      1.000        0.000        0.001        0.000
## 17       D3      0.998      1.000        0.000        0.001        0.000
## 18      D33      0.998      1.000        0.000        0.001        0.000
## 19      D34      0.999      1.000        0.000        0.001        0.000
## 20      D35      0.998      1.000        0.000        0.001        0.000
## 21      D37      0.999      1.000        0.000        0.001        0.000
## 22      D38      0.998      1.000        0.000        0.001        0.000
## 23      D39      0.998      1.000        0.000        0.001        0.000
## 24       D4      0.998      1.000        0.000        0.001        0.000
## 25      D40      0.998      1.000        0.000        0.001        0.000
## 26      D41      0.998      1.000        0.000        0.001        0.000
## 27      D42      0.998      1.000        0.000        0.001        0.000
## 28      D43      0.998      1.000        0.000        0.001        0.000
## 29      D44      0.999      1.000        0.000        0.001        0.000
## 30      D46      0.997      1.000        0.000        0.001        0.000
## 31      D47      0.998      1.000        0.000        0.001        0.000
## 32      D48      0.999      1.000        0.000        0.001        0.000
## 33      D49      0.998      1.000        0.000        0.001        0.000
## 34       D5      0.999      1.000        0.000        0.001        0.000
## 35      D50      0.999      1.000        0.000        0.001        0.000
## 36      D51      0.998      1.000        0.000        0.001        0.000
## 37      D52      0.999      1.000        0.000        0.001        0.000
## 38      D53      0.999      1.000        0.000        0.001        0.000
## 39      D54      0.998      1.000        0.000        0.001        0.000
## 40      D55      0.998      1.000        0.000        0.001        0.000
## 41      D56      0.999      1.000        0.000        0.001        0.000
## 42      D57      0.998      1.000        0.000        0.001        0.000
## 43      D58      0.998      1.000        0.000        0.001        0.000
## 44      D59      0.998      1.000        0.000        0.001        0.000
## 45       D6      0.999      1.000        0.000        0.001        0.000
## 46      D60      0.998      1.000        0.000        0.001        0.000
## 47      D61      0.999      1.000        0.000        0.001        0.000
## 48      D62      0.998      1.000        0.000        0.001        0.000
## 49      D63      0.998      1.000        0.000        0.001        0.000
## 50      D64      0.998      1.000        0.000        0.001        0.000
## 51  D65_rep      0.970      1.000        0.000        0.012        0.000
## 52       D7      0.999      1.000        0.000        0.001        0.000
## 53  D71_rep      0.988      1.000        0.000        0.001        0.000
## 54      D72      0.997      1.000        0.000        0.001        0.000
## 55      D73      0.998      1.000        0.000        0.001        0.000
## 56      D75      0.998      1.000        0.000        0.001        0.000
## 57      D76      0.998      1.000        0.000        0.001        0.000
## 58      D77      0.998      1.000        0.000        0.001        0.000
## 59      D78      0.999      1.000        0.000        0.001        0.000
## 60      D79      0.998      1.000        0.000        0.001        0.000
## 61       D8      0.999      1.000        0.000        0.001        0.000
## 62      D80      0.999      1.000        0.000        0.001        0.000
## 63      D81      0.999      1.000        0.000        0.001        0.000
## 64      D82      0.998      1.000        0.000        0.001        0.000
## 65      D83      0.999      1.000        0.000        0.001        0.000
## 66      D84      0.998      1.000        0.000        0.001        0.000
## 67      D85      0.998      1.000        0.000        0.001        0.000
## 68      D86      0.998      1.000        0.000        0.001        0.000
## 69      D87      0.998      1.000        0.000        0.001        0.000
## 70      D88      0.998      1.000        0.000        0.001        0.000
## 71      D89      0.998      1.000        0.000        0.001        0.000
## 72      D90      0.997      1.000        0.000        0.001        0.000
## 73      D91      0.998      1.000        0.000        0.001        0.000
## 74      D92      0.997      1.000        0.000        0.001        0.000
## 75      D93      0.998      1.000        0.000        0.001        0.000
## 76      D94      0.998      1.000        0.000        0.001        0.000
## 77      J14      0.000      0.001        0.998        1.000        0.000
## 78      J15      0.000      0.001        0.997        1.000        0.000
## 79      J16      0.000      0.001        0.998        1.000        0.000
## 80      J17      0.000      0.001        0.999        1.000        0.000
## 81      J18      0.000      0.024        0.975        1.000        0.000
## 82      J19      0.000      0.001        0.999        1.000        0.000
## 83      J20      0.000      0.001        0.999        1.000        0.000
## 84      J21      0.000      0.001        0.999        1.000        0.000
## 85      J22      0.000      0.001        0.999        1.000        0.000
## 86      J23      0.000      0.001        0.998        1.000        0.000
## 87      J24      0.000      0.001        0.999        1.000        0.000
## 88      J25      0.000      0.001        0.999        1.000        0.000
## 89      J26      0.000      0.001        0.999        1.000        0.000
## 90      J27      0.000      0.001        0.999        1.000        0.000
## 91       K1      0.000      0.001        0.999        1.000        0.000
## 92       K2      0.000      0.001        0.999        1.000        0.000
## 93       K3      0.000      0.001        0.999        1.000        0.000
## 94      K32      0.000      0.001        0.999        1.000        0.000
## 95      K33      0.000      0.001        0.999        1.000        0.000
## 96      K34      0.000      0.001        0.998        1.000        0.000
## 97       K4      0.000      0.001        0.999        1.000        0.000
## 98       K5      0.000      0.001        0.998        1.000        0.000
## 99       L1      0.000      0.001        0.998        1.000        0.000
## 100     L10      0.000      0.001        0.991        1.000        0.000
## 101     L11      0.000      0.001        0.999        1.000        0.000
## 102     L12      0.000      0.001        0.998        1.000        0.000
## 103     L13      0.000      0.001        0.998        1.000        0.000
## 104     L14      0.000      0.001        0.999        1.000        0.000
## 105     L15      0.000      0.001        0.998        1.000        0.000
## 106     L16      0.000      0.001        0.999        1.000        0.000
## 107     L17      0.000      0.001        0.999        1.000        0.000
## 108 L18_rep      0.000      0.001        0.998        1.000        0.000
## 109     L19      0.000      0.001        0.999        1.000        0.000
## 110      L2      0.000      0.001        0.998        1.000        0.000
## 111     L20      0.000      0.001        0.999        1.000        0.000
## 112     L21      0.000      0.001        0.998        1.000        0.000
## 113     L22      0.000      0.001        0.999        1.000        0.000
## 114     L23      0.000      0.001        0.999        1.000        0.000
## 115     L24      0.000      0.001        0.998        1.000        0.000
## 116     L26      0.000      0.001        0.998        1.000        0.000
## 117     L27      0.000      0.001        0.998        1.000        0.000
## 118     L28      0.000      0.001        0.999        1.000        0.000
## 119      L3      0.000      0.001        0.998        1.000        0.000
## 120      L4      0.000      0.001        0.999        1.000        0.000
## 121      L5      0.000      0.001        0.998        1.000        0.000
## 122      L6      0.000      0.001        0.999        1.000        0.000
## 123      L7      0.000      0.001        0.995        1.000        0.000
## 124      L8      0.000      0.001        0.999        1.000        0.000
## 125  L9_rep      0.000      0.001        0.998        1.000        0.000
## 126      N1      0.000      0.001        0.995        1.000        0.000
## 127     N10      0.000      0.002        0.981        1.000        0.000
## 128     N11      0.000      0.001        0.997        1.000        0.000
## 129     N12      0.000      0.001        0.999        1.000        0.000
## 130     N13      0.000      0.001        0.999        1.000        0.000
## 131     N14      0.000      0.001        0.999        1.000        0.000
## 132     N15      0.000      0.001        0.999        1.000        0.000
## 133     N16      0.000      0.001        0.998        1.000        0.000
## 134     N17      0.000      0.001        0.999        1.000        0.000
## 135     N18      0.000      0.002        0.973        1.000        0.000
## 136     N19      0.000      0.001        0.998        1.000        0.000
## 137     N20      0.000      0.001        0.999        1.000        0.000
## 138     N21      0.000      0.001        0.998        1.000        0.000
## 139     N24      0.000      0.001        0.998        1.000        0.000
## 140     N25      0.000      0.001        0.998        1.000        0.000
## 141     N26      0.000      0.001        0.972        1.000        0.000
## 142     N27      0.000      0.001        0.987        1.000        0.000
## 143     N28      0.000      0.001        0.999        1.000        0.000
## 144      N3      0.000      0.002        0.987        1.000        0.000
## 145      N4      0.000      0.001        0.999        1.000        0.000
## 146      N5      0.000      0.001        0.998        1.000        0.000
## 147      N6      0.000      0.001        0.999        1.000        0.000
## 148      N8      0.000      0.001        0.998        1.000        0.000
## 149      N9      0.000      0.001        0.998        1.000        0.000
## 150     R10      0.000      0.001        0.999        1.000        0.000
## 151     R11      0.000      0.001        0.996        1.000        0.000
## 152     R12      0.000      0.001        0.998        1.000        0.000
## 153     R13      0.000      0.001        0.997        1.000        0.000
## 154     R14      0.000      0.001        0.999        1.000        0.000
## 155     R15      0.000      0.001        0.998        1.000        0.000
## 156     R16      0.000      0.001        0.999        1.000        0.000
## 157     R17      0.000      0.001        0.999        1.000        0.000
## 158     R18      0.000      0.001        0.999        1.000        0.000
## 159     R19      0.000      0.001        0.979        1.000        0.000
## 160      R2      0.000      0.001        0.998        1.000        0.000
## 161     R20      0.000      0.001        0.998        1.000        0.000
## 162      R3      0.000      0.001        0.998        1.000        0.000
## 163      R5      0.000      0.001        0.999        1.000        0.000
## 164      R6      0.000      0.001        0.999        1.000        0.000
## 165      R7      0.000      0.001        0.999        1.000        0.000
## 166      R9      0.000      0.001        0.999        1.000        0.000
## 167      S1      0.000      0.001        0.999        1.000        0.000
## 168 S11_rep      0.000      0.001        0.998        1.000        0.000
## 169     S12      0.000      0.001        0.998        1.000        0.000
## 170     S13      0.000      0.001        0.998        1.000        0.000
## 171     S15      0.000      0.001        0.999        1.000        0.000
## 172     S16      0.000      0.001        0.999        1.000        0.000
## 173     S17      0.000      0.001        0.999        1.000        0.000
## 174     S18      0.000      0.001        0.998        1.000        0.000
## 175     S19      0.000      0.001        0.999        1.000        0.000
## 176      S2      0.000      0.001        0.999        1.000        0.000
## 177     S20      0.000      0.001        0.998        1.000        0.000
## 178     S21      0.000      0.001        0.999        1.000        0.000
## 179     S22      0.000      0.001        0.998        1.000        0.000
## 180     S23      0.000      0.001        0.999        1.000        0.000
## 181     S24      0.000      0.001        0.999        1.000        0.000
## 182      S3      0.000      0.001        0.993        1.000        0.000
## 183      S4      0.000      0.001        0.999        1.000        0.000
## 184      S6      0.000      0.001        0.999        1.000        0.000
## 185      S7      0.000      0.001        0.998        1.000        0.000
## 186      S8      0.000      0.001        0.999        1.000        0.000
## 187      S9      0.000      0.001        0.998        1.000        0.000
##     carinulata_H carinata_L carinata_H elongata sublineata carinulata carinata
## 1          0.112          0      0.001    0.000      0.912      0.088    0.000
## 2          0.089          0      0.001    0.931      0.000      0.068    0.000
## 3          0.046          0      0.001    0.970      0.000      0.030    0.000
## 4          0.001          0      0.001    1.000      0.000      0.000    0.000
## 5          0.001          0      0.001    1.000      0.000      0.000    0.000
## 6          0.001          0      0.001    1.000      0.000      0.000    0.000
## 7          0.001          0      0.001    1.000      0.000      0.000    0.000
## 8          0.001          0      0.005    0.998      0.001      0.000    0.001
## 9          0.001          0      0.001    1.000      0.000      0.000    0.000
## 10         0.001          0      0.001    1.000      0.000      0.000    0.000
## 11         0.001          0      0.001    1.000      0.000      0.000    0.000
## 12         0.001          0      0.001    1.000      0.000      0.000    0.000
## 13         0.001          0      0.001    1.000      0.000      0.000    0.000
## 14         0.001          0      0.001    1.000      0.000      0.000    0.000
## 15         0.001          0      0.001    1.000      0.000      0.000    0.000
## 16         0.001          0      0.001    1.000      0.000      0.000    0.000
## 17         0.001          0      0.001    1.000      0.000      0.000    0.000
## 18         0.001          0      0.001    1.000      0.000      0.000    0.000
## 19         0.001          0      0.001    1.000      0.000      0.000    0.000
## 20         0.001          0      0.001    1.000      0.000      0.000    0.000
## 21         0.001          0      0.001    1.000      0.000      0.000    0.000
## 22         0.001          0      0.001    1.000      0.000      0.000    0.000
## 23         0.001          0      0.001    1.000      0.000      0.000    0.000
## 24         0.001          0      0.001    1.000      0.000      0.000    0.000
## 25         0.001          0      0.001    1.000      0.000      0.000    0.000
## 26         0.001          0      0.001    1.000      0.000      0.000    0.000
## 27         0.001          0      0.001    1.000      0.000      0.000    0.000
## 28         0.001          0      0.001    1.000      0.000      0.000    0.000
## 29         0.001          0      0.001    1.000      0.000      0.000    0.000
## 30         0.001          0      0.001    1.000      0.000      0.000    0.000
## 31         0.001          0      0.001    1.000      0.000      0.000    0.000
## 32         0.001          0      0.001    1.000      0.000      0.000    0.000
## 33         0.001          0      0.001    1.000      0.000      0.000    0.000
## 34         0.001          0      0.001    1.000      0.000      0.000    0.000
## 35         0.001          0      0.001    1.000      0.000      0.000    0.000
## 36         0.001          0      0.001    1.000      0.000      0.000    0.000
## 37         0.001          0      0.001    1.000      0.000      0.000    0.000
## 38         0.001          0      0.001    1.000      0.000      0.000    0.000
## 39         0.001          0      0.001    1.000      0.000      0.000    0.000
## 40         0.001          0      0.001    1.000      0.000      0.000    0.000
## 41         0.001          0      0.001    1.000      0.000      0.000    0.000
## 42         0.001          0      0.001    1.000      0.000      0.000    0.000
## 43         0.001          0      0.001    1.000      0.000      0.000    0.000
## 44         0.001          0      0.001    1.000      0.000      0.000    0.000
## 45         0.001          0      0.001    1.000      0.000      0.000    0.000
## 46         0.001          0      0.001    1.000      0.000      0.000    0.000
## 47         0.001          0      0.001    1.000      0.000      0.000    0.000
## 48         0.001          0      0.001    1.000      0.000      0.000    0.000
## 49         0.001          0      0.001    1.000      0.000      0.000    0.000
## 50         0.001          0      0.001    1.000      0.000      0.000    0.000
## 51         0.029          0      0.001    0.989      0.002      0.010    0.000
## 52         0.001          0      0.001    1.000      0.000      0.000    0.000
## 53         0.012          0      0.001    0.998      0.000      0.002    0.000
## 54         0.001          0      0.001    1.000      0.000      0.000    0.000
## 55         0.001          0      0.001    1.000      0.000      0.000    0.000
## 56         0.001          0      0.001    1.000      0.000      0.000    0.000
## 57         0.001          0      0.001    1.000      0.000      0.000    0.000
## 58         0.001          0      0.001    1.000      0.000      0.000    0.000
## 59         0.001          0      0.001    1.000      0.000      0.000    0.000
## 60         0.001          0      0.001    1.000      0.000      0.000    0.000
## 61         0.001          0      0.001    1.000      0.000      0.000    0.000
## 62         0.001          0      0.001    1.000      0.000      0.000    0.000
## 63         0.001          0      0.001    1.000      0.000      0.000    0.000
## 64         0.001          0      0.001    1.000      0.000      0.000    0.000
## 65         0.001          0      0.001    1.000      0.000      0.000    0.000
## 66         0.001          0      0.001    1.000      0.000      0.000    0.000
## 67         0.001          0      0.001    1.000      0.000      0.000    0.000
## 68         0.001          0      0.001    1.000      0.000      0.000    0.000
## 69         0.001          0      0.001    1.000      0.000      0.000    0.000
## 70         0.001          0      0.001    1.000      0.000      0.000    0.000
## 71         0.001          0      0.001    1.000      0.000      0.000    0.000
## 72         0.001          0      0.001    1.000      0.000      0.000    0.000
## 73         0.001          0      0.001    1.000      0.000      0.000    0.000
## 74         0.001          0      0.001    1.000      0.000      0.000    0.000
## 75         0.001          0      0.001    1.000      0.000      0.000    0.000
## 76         0.001          0      0.001    1.000      0.000      0.000    0.000
## 77         0.001          0      0.001    0.000      1.000      0.000    0.000
## 78         0.001          0      0.002    0.000      1.000      0.000    0.000
## 79         0.001          0      0.001    0.000      1.000      0.000    0.000
## 80         0.001          0      0.001    0.000      1.000      0.000    0.000
## 81         0.001          0      0.001    0.012      0.988      0.000    0.000
## 82         0.001          0      0.001    0.000      1.000      0.000    0.000
## 83         0.001          0      0.001    0.000      1.000      0.000    0.000
## 84         0.001          0      0.001    0.000      1.000      0.000    0.000
## 85         0.001          0      0.001    0.000      1.000      0.000    0.000
## 86         0.001          0      0.001    0.000      1.000      0.000    0.000
## 87         0.001          0      0.001    0.000      1.000      0.000    0.000
## 88         0.001          0      0.001    0.000      1.000      0.000    0.000
## 89         0.001          0      0.001    0.000      1.000      0.000    0.000
## 90         0.001          0      0.001    0.000      1.000      0.000    0.000
## 91         0.001          0      0.001    0.000      1.000      0.000    0.000
## 92         0.001          0      0.001    0.000      1.000      0.000    0.000
## 93         0.001          0      0.001    0.000      1.000      0.000    0.000
## 94         0.001          0      0.001    0.000      1.000      0.000    0.000
## 95         0.001          0      0.001    0.000      1.000      0.000    0.000
## 96         0.001          0      0.001    0.000      1.000      0.000    0.000
## 97         0.001          0      0.001    0.000      1.000      0.000    0.000
## 98         0.001          0      0.001    0.000      1.000      0.000    0.000
## 99         0.001          0      0.001    0.000      1.000      0.000    0.000
## 100        0.004          0      0.005    0.000      0.999      0.000    0.001
## 101        0.001          0      0.001    0.000      1.000      0.000    0.000
## 102        0.001          0      0.001    0.000      1.000      0.000    0.000
## 103        0.001          0      0.001    0.000      1.000      0.000    0.000
## 104        0.001          0      0.001    0.000      1.000      0.000    0.000
## 105        0.001          0      0.001    0.000      1.000      0.000    0.000
## 106        0.001          0      0.001    0.000      1.000      0.000    0.000
## 107        0.001          0      0.001    0.000      1.000      0.000    0.000
## 108        0.001          0      0.001    0.000      1.000      0.000    0.000
## 109        0.001          0      0.001    0.000      1.000      0.000    0.000
## 110        0.001          0      0.001    0.000      1.000      0.000    0.000
## 111        0.001          0      0.001    0.000      1.000      0.000    0.000
## 112        0.001          0      0.001    0.000      1.000      0.000    0.000
## 113        0.001          0      0.001    0.000      1.000      0.000    0.000
## 114        0.001          0      0.001    0.000      1.000      0.000    0.000
## 115        0.001          0      0.001    0.000      1.000      0.000    0.000
## 116        0.001          0      0.001    0.000      1.000      0.000    0.000
## 117        0.001          0      0.001    0.000      1.000      0.000    0.000
## 118        0.001          0      0.001    0.000      1.000      0.000    0.000
## 119        0.001          0      0.001    0.000      1.000      0.000    0.000
## 120        0.001          0      0.001    0.000      1.000      0.000    0.000
## 121        0.001          0      0.001    0.000      1.000      0.000    0.000
## 122        0.001          0      0.001    0.000      1.000      0.000    0.000
## 123        0.001          0      0.005    0.000      0.999      0.000    0.001
## 124        0.001          0      0.001    0.000      1.000      0.000    0.000
## 125        0.001          0      0.001    0.000      1.000      0.000    0.000
## 126        0.001          0      0.004    0.000      0.999      0.000    0.001
## 127        0.018          0      0.002    0.000      0.994      0.006    0.000
## 128        0.002          0      0.001    0.000      1.000      0.000    0.000
## 129        0.001          0      0.001    0.000      1.000      0.000    0.000
## 130        0.001          0      0.001    0.000      1.000      0.000    0.000
## 131        0.001          0      0.001    0.000      1.000      0.000    0.000
## 132        0.001          0      0.001    0.000      1.000      0.000    0.000
## 133        0.001          0      0.001    0.000      1.000      0.000    0.000
## 134        0.001          0      0.001    0.000      1.000      0.000    0.000
## 135        0.023          0      0.017    0.000      0.989      0.008    0.002
## 136        0.001          0      0.001    0.000      1.000      0.000    0.000
## 137        0.001          0      0.001    0.000      1.000      0.000    0.000
## 138        0.001          0      0.001    0.000      1.000      0.000    0.000
## 139        0.001          0      0.001    0.000      1.000      0.000    0.000
## 140        0.001          0      0.001    0.000      1.000      0.000    0.000
## 141        0.027          0      0.004    0.000      0.986      0.013    0.001
## 142        0.013          0      0.001    0.000      0.997      0.002    0.000
## 143        0.001          0      0.001    0.000      1.000      0.000    0.000
## 144        0.012          0      0.001    0.000      0.998      0.002    0.000
## 145        0.001          0      0.001    0.000      1.000      0.000    0.000
## 146        0.001          0      0.001    0.000      1.000      0.000    0.000
## 147        0.001          0      0.001    0.000      1.000      0.000    0.000
## 148        0.001          0      0.001    0.000      1.000      0.000    0.000
## 149        0.001          0      0.001    0.000      1.000      0.000    0.000
## 150        0.001          0      0.001    0.000      1.000      0.000    0.000
## 151        0.001          0      0.003    0.000      0.999      0.000    0.000
## 152        0.001          0      0.001    0.000      1.000      0.000    0.000
## 153        0.001          0      0.002    0.000      0.999      0.000    0.000
## 154        0.001          0      0.001    0.000      1.000      0.000    0.000
## 155        0.001          0      0.001    0.000      1.000      0.000    0.000
## 156        0.001          0      0.001    0.000      1.000      0.000    0.000
## 157        0.001          0      0.001    0.000      1.000      0.000    0.000
## 158        0.001          0      0.001    0.000      1.000      0.000    0.000
## 159        0.016          0      0.013    0.000      0.994      0.004    0.002
## 160        0.001          0      0.001    0.000      1.000      0.000    0.000
## 161        0.001          0      0.001    0.000      1.000      0.000    0.000
## 162        0.001          0      0.001    0.000      1.000      0.000    0.000
## 163        0.001          0      0.001    0.000      1.000      0.000    0.000
## 164        0.001          0      0.001    0.000      1.000      0.000    0.000
## 165        0.001          0      0.001    0.000      1.000      0.000    0.000
## 166        0.001          0      0.001    0.000      1.000      0.000    0.000
## 167        0.001          0      0.001    0.000      1.000      0.000    0.000
## 168        0.001          0      0.001    0.000      1.000      0.000    0.000
## 169        0.001          0      0.001    0.000      1.000      0.000    0.000
## 170        0.001          0      0.001    0.000      1.000      0.000    0.000
## 171        0.001          0      0.001    0.000      1.000      0.000    0.000
## 172        0.001          0      0.001    0.000      1.000      0.000    0.000
## 173        0.001          0      0.001    0.000      1.000      0.000    0.000
## 174        0.001          0      0.001    0.000      1.000      0.000    0.000
## 175        0.001          0      0.001    0.000      1.000      0.000    0.000
## 176        0.001          0      0.001    0.000      1.000      0.000    0.000
## 177        0.001          0      0.001    0.000      1.000      0.000    0.000
## 178        0.001          0      0.001    0.000      1.000      0.000    0.000
## 179        0.001          0      0.001    0.000      1.000      0.000    0.000
## 180        0.001          0      0.001    0.000      1.000      0.000    0.000
## 181        0.001          0      0.001    0.000      1.000      0.000    0.000
## 182        0.002          0      0.003    0.000      0.999      0.000    0.001
## 183        0.001          0      0.001    0.000      1.000      0.000    0.000
## 184        0.001          0      0.001    0.000      1.000      0.000    0.000
## 185        0.001          0      0.001    0.000      1.000      0.000    0.000
## 186        0.001          0      0.001    0.000      1.000      0.000    0.000
## 187        0.001          0      0.001    0.000      1.000      0.000    0.000
##       popID
## 1      47TX
## 2      44GR
## 3      41CR
## 4      38CR
## 5      41CR
## 6      41CR
## 7      41CR
## 8      43GR
## 9      43GR
## 10     43GR
## 11     38CR
## 12     43GR
## 13     44GR
## 14     44GR
## 15     44GR
## 16     44GR
## 17     38CR
## 18     43GR
## 19     43GR
## 20     43GR
## 21     44GR
## 22     44GR
## 23     44GR
## 24     38CR
## 25     44GR
## 26     41CR
## 27     41CR
## 28     41CR
## 29     41CR
## 30     37CR
## 31     37CR
## 32     37CR
## 33     38CR
## 34     39CR
## 35     38CR
## 36     38CR
## 37     38CR
## 38     39CR
## 39     39CR
## 40     39CR
## 41     39CR
## 42     41CR
## 43     41CR
## 44     41CR
## 45     39CR
## 46     41CR
## 47     43GR
## 48     43GR
## 49     43GR
## 50     43GR
## 51     44GR
## 52     39CR
## 53     37CR
## 54     37CR
## 55     38CR
## 56     38CR
## 57     38CR
## 58     39CR
## 59     39CR
## 60     39CR
## 61     39CR
## 62     39CR
## 63     41CR
## 64     41CR
## 65     41CR
## 66     41CR
## 67     43GR
## 68     43GR
## 69     43GR
## 70     43GR
## 71     44GR
## 72     44GR
## 73     44GR
## 74     44GR
## 75     44GR
## 76     44GR
## 77     51TX
## 78     51TX
## 79     51TX
## 80     51TX
## 81     51TX
## 82     51TX
## 83     51TX
## 84     51TX
## 85     51TX
## 86     51TX
## 87     51TX
## 88     51TX
## 89     51TX
## 90     51TX
## 91     50TX
## 92     50TX
## 93     50TX
## 94     50TX
## 95     50TX
## 96     50TX
## 97     50TX
## 98     50TX
## 99     49TX
## 100    49TX
## 101    49TX
## 102    49TX
## 103    49TX
## 104    49TX
## 105    49TX
## 106    49TX
## 107    49TX
## 108    49TX
## 109    49TX
## 110    49TX
## 111    49TX
## 112    49TX
## 113    49TX
## 114    49TX
## 115    49TX
## 116    49TX
## 117    49TX
## 118    49TX
## 119    49TX
## 120    49TX
## 121    49TX
## 122    49TX
## 123    49TX
## 124    49TX
## 125    49TX
## 126    47TX
## 127    47TX
## 128    47TX
## 129    47TX
## 130    47TX
## 131    47TX
## 132    47TX
## 133    47TX
## 134    47TX
## 135    47TX
## 136    47TX
## 137    47TX
## 138    47TX
## 139    47TX
## 140    47TX
## 141    47TX
## 142    47TX
## 143    47TX
## 144    47TX
## 145    47TX
## 146    47TX
## 147    47TX
## 148    47TX
## 149    47TX
## 150 SUB_LAB
## 151 SUB_LAB
## 152 SUB_LAB
## 153 SUB_LAB
## 154 SUB_LAB
## 155 SUB_LAB
## 156 SUB_LAB
## 157 SUB_LAB
## 158 SUB_LAB
## 159 SUB_LAB
## 160 SUB_LAB
## 161 SUB_LAB
## 162 SUB_LAB
## 163 SUB_LAB
## 164 SUB_LAB
## 165 SUB_LAB
## 166 SUB_LAB
## 167    52TX
## 168    52TX
## 169    52TX
## 170    52TX
## 171    52TX
## 172    52TX
## 173    52TX
## 174    52TX
## 175    52TX
## 176    52TX
## 177    52TX
## 178    52TX
## 179    52TX
## 180    52TX
## 181    52TX
## 182    52TX
## 183    52TX
## 184    52TX
## 185    52TX
## 186    52TX
## 187    52TX
find_threshold <- ci %>%
  semi_join(diagnostic) %>%
  arrange(desc(carinulata_L))
## Joining, by = "popID"
t <- find_threshold[1,"carinulata"]
t
## [1] 0.088
t_low <- find_threshold[1,"carinulata_L"]
t_low
## [1] 0.066

Each hybrid type

How many hybrids with carinulata?

carinu_hyb <- ci %>%
  dplyr::filter(carinulata_L > t_low ) %>%
  dplyr::filter(elongata_L > t_low |
           sublineata_L > t_low | 
             carinata_L > t_low) %>%
  arrange(desc(carinulata_L))

carinu_hyb
##    sampID elongata_L elongata_H sublineata_L sublineata_H carinulata_L
## 1  G5_rep          0      0.001        0.031        0.116        0.633
## 2 G48_rep          0      0.001        0.890        0.932        0.068
##   carinulata_H carinata_L carinata_H elongata sublineata carinulata carinata
## 1        0.695      0.216      0.312        0      0.071      0.664    0.265
## 2        0.109      0.000      0.001        0      0.912      0.088    0.000
##   popID
## 1  28NM
## 2   6KS

How many tri-specific hybrids?

ci %>%
  dplyr::filter(carinata_L > t_low & elongata_L > t_low & sublineata_L > t_low) 
##  [1] sampID       elongata_L   elongata_H   sublineata_L sublineata_H
##  [6] carinulata_L carinulata_H carinata_L   carinata_H   elongata    
## [11] sublineata   carinulata   carinata     popID       
## <0 rows> (or 0-length row.names)
## If we count less-conservatively:
ci %>%
  dplyr::filter(carinata_L >.01 & elongata_L> .01 & sublineata_L > .01) %>%
  arrange(desc(elongata_L))
##   sampID elongata_L elongata_H sublineata_L sublineata_H carinulata_L
## 1    H10      0.024      0.079        0.140        0.258        0.015
## 2    F84      0.018      0.052        0.867        0.940        0.000
## 3     G2      0.018      0.051        0.111        0.202        0.000
##   carinulata_H carinata_L carinata_H elongata sublineata carinulata carinata
## 1        0.056      0.655      0.777    0.050      0.198      0.035    0.717
## 2        0.001      0.024      0.100    0.034      0.905      0.000    0.060
## 3        0.001      0.762      0.856    0.034      0.155      0.000    0.811
##   popID
## 1  28NM
## 2  26NM
## 3  15TX

How many hybrids with sublineata?

ci %>%
  dplyr::filter(sublineata_L > t_low) %>%
  dplyr::filter(carinulata_L > t_low |
           elongata_L > t_low | carinata_L > t_low) %>%
  arrange(desc(sublineata_L)) 
##     sampID elongata_L elongata_H sublineata_L sublineata_H carinulata_L
## 1  G48_rep      0.000      0.001        0.890        0.932        0.068
## 2      G10      0.088      0.133        0.866        0.912        0.000
## 3      G50      0.092      0.140        0.860        0.908        0.000
## 4      F95      0.000      0.001        0.844        0.928        0.000
## 5      G11      0.000      0.001        0.818        0.912        0.000
## 6       G6      0.000      0.001        0.815        0.897        0.000
## 7      H21      0.000      0.001        0.797        0.887        0.000
## 8      H71      0.000      0.036        0.796        0.885        0.000
## 9      H23      0.000      0.037        0.793        0.878        0.000
## 10     G63      0.000      0.019        0.762        0.869        0.000
## 11     H26      0.198      0.263        0.736        0.802        0.000
## 12     F94      0.002      0.051        0.727        0.838        0.000
## 13     H12      0.000      0.001        0.727        0.828        0.000
## 14      G9      0.000      0.001        0.721        0.820        0.000
## 15     G52      0.000      0.021        0.667        0.771        0.000
## 16     G83      0.000      0.002        0.642        0.802        0.000
## 17     F96      0.000      0.003        0.619        0.739        0.000
## 18     G19      0.000      0.001        0.572        0.704        0.000
## 19     H67      0.000      0.001        0.472        0.599        0.000
## 20     H68      0.000      0.001        0.397        0.521        0.000
## 21     F92      0.000      0.001        0.380        0.511        0.000
## 22     H22      0.000      0.007        0.378        0.503        0.000
## 23      H9      0.000      0.002        0.348        0.466        0.000
## 24     H64      0.000      0.001        0.291        0.412        0.000
## 25     G38      0.000      0.001        0.271        0.393        0.000
## 26     G78      0.000      0.001        0.154        0.262        0.000
## 27     H10      0.024      0.079        0.140        0.258        0.015
## 28      G2      0.018      0.051        0.111        0.202        0.000
##    carinulata_H carinata_L carinata_H elongata sublineata carinulata carinata
## 1         0.109      0.000      0.001    0.000      0.912      0.088    0.000
## 2         0.001      0.000      0.003    0.110      0.890      0.000    0.001
## 3         0.001      0.000      0.001    0.115      0.885      0.000    0.000
## 4         0.001      0.071      0.156    0.000      0.888      0.000    0.111
## 5         0.001      0.088      0.182    0.000      0.867      0.000    0.133
## 6         0.001      0.103      0.185    0.000      0.858      0.000    0.142
## 7         0.001      0.113      0.202    0.000      0.844      0.000    0.156
## 8         0.001      0.096      0.190    0.017      0.843      0.000    0.140
## 9         0.001      0.100      0.189    0.020      0.837      0.000    0.143
## 10        0.004      0.125      0.235    0.003      0.817      0.000    0.179
## 11        0.001      0.000      0.001    0.230      0.770      0.000    0.000
## 12        0.001      0.130      0.247    0.028      0.785      0.000    0.187
## 13        0.001      0.172      0.272    0.000      0.779      0.000    0.221
## 14        0.001      0.180      0.279    0.000      0.772      0.000    0.228
## 15        0.001      0.223      0.330    0.004      0.720      0.000    0.276
## 16        0.001      0.198      0.357    0.000      0.725      0.000    0.274
## 17        0.001      0.260      0.381    0.001      0.680      0.000    0.320
## 18        0.001      0.296      0.428    0.000      0.639      0.000    0.360
## 19        0.001      0.401      0.527    0.000      0.536      0.000    0.464
## 20        0.001      0.479      0.603    0.000      0.459      0.000    0.541
## 21        0.001      0.489      0.620    0.000      0.445      0.000    0.555
## 22        0.001      0.496      0.621    0.001      0.440      0.000    0.559
## 23        0.001      0.533      0.651    0.000      0.407      0.000    0.593
## 24        0.001      0.588      0.709    0.000      0.350      0.000    0.650
## 25        0.001      0.607      0.729    0.000      0.331      0.000    0.668
## 26        0.001      0.737      0.846    0.000      0.207      0.000    0.793
## 27        0.056      0.655      0.777    0.050      0.198      0.035    0.717
## 28        0.001      0.762      0.856    0.034      0.155      0.000    0.811
##    popID
## 1    6KS
## 2   26NM
## 3   20NM
## 4   26NM
## 5   26NM
## 6   28NM
## 7   26NM
## 8   26NM
## 9   26NM
## 10  20NM
## 11  21NM
## 12  26NM
## 13  26NM
## 14  26NM
## 15  20NM
## 16  21NM
## 17  26NM
## 18  27NM
## 19  28NM
## 20  28NM
## 21  28NM
## 22  26NM
## 23  28NM
## 24  18TX
## 25  20NM
## 26  16TX
## 27  28NM
## 28  15TX
ci %>%
  dplyr::filter(sublineata_L > t) %>%
  dplyr::filter(carinulata_L > t_low |
           elongata_L > t_low | carinata_L > t_low) %>%
  arrange(desc(sublineata_L)) %>%
  count()
##    n
## 1 28

D. elongata hybrids

ci %>%
  dplyr::filter(elongata_L > t_low) %>%
  dplyr::filter(carinulata_L > t_low |
           sublineata_L > t_low | carinata_L > t_low) %>%
  arrange(desc(elongata_L))
##    sampID elongata_L elongata_H sublineata_L sublineata_H carinulata_L
## 1     E29      0.221      0.305        0.000        0.002            0
## 2     H26      0.198      0.263        0.736        0.802            0
## 3     G50      0.092      0.140        0.860        0.908            0
## 4     G10      0.088      0.133        0.866        0.912            0
## 5     G87      0.082      0.137        0.000        0.005            0
## 6 G01_rep      0.078      0.129        0.000        0.001            0
##   carinulata_H carinata_L carinata_H elongata sublineata carinulata carinata
## 1        0.001      0.694      0.779    0.262      0.000          0    0.737
## 2        0.001      0.000      0.001    0.230      0.770          0    0.000
## 3        0.001      0.000      0.001    0.115      0.885          0    0.000
## 4        0.001      0.000      0.003    0.110      0.890          0    0.001
## 5        0.001      0.862      0.917    0.109      0.001          0    0.890
## 6        0.001      0.870      0.922    0.103      0.000          0    0.897
##   popID
## 1  15TX
## 2  21NM
## 3  20NM
## 4  26NM
## 5  15TX
## 6  15TX
ci %>%
  dplyr::filter(elongata_L > t_low) %>%
  dplyr::filter(carinulata_L > t_low |
           sublineata_L > t_low | carinata_L > t_low) %>%
  arrange(desc(elongata_L)) %>%
  count()
##   n
## 1 6

Specific bi-parental types

How many carinata-sublineata hybrids? How many carinata-elongata hybrids? How many elong-sublineata hybrids?

car_sub_hyb <- ci %>%
  dplyr::filter(carinata_L >t_low  & sublineata_L > t_low & elongata_L < t_low)
car_sub_hyb
##    sampID elongata_L elongata_H sublineata_L sublineata_H carinulata_L
## 1     F92      0.000      0.001        0.380        0.511        0.000
## 2     F94      0.002      0.051        0.727        0.838        0.000
## 3     F95      0.000      0.001        0.844        0.928        0.000
## 4     F96      0.000      0.003        0.619        0.739        0.000
## 5     G11      0.000      0.001        0.818        0.912        0.000
## 6     G19      0.000      0.001        0.572        0.704        0.000
## 7      G2      0.018      0.051        0.111        0.202        0.000
## 8     G38      0.000      0.001        0.271        0.393        0.000
## 9     G52      0.000      0.021        0.667        0.771        0.000
## 10     G6      0.000      0.001        0.815        0.897        0.000
## 11    G63      0.000      0.019        0.762        0.869        0.000
## 12    G78      0.000      0.001        0.154        0.262        0.000
## 13    G83      0.000      0.002        0.642        0.802        0.000
## 14     G9      0.000      0.001        0.721        0.820        0.000
## 15    H10      0.024      0.079        0.140        0.258        0.015
## 16    H12      0.000      0.001        0.727        0.828        0.000
## 17    H21      0.000      0.001        0.797        0.887        0.000
## 18    H22      0.000      0.007        0.378        0.503        0.000
## 19    H23      0.000      0.037        0.793        0.878        0.000
## 20    H64      0.000      0.001        0.291        0.412        0.000
## 21    H67      0.000      0.001        0.472        0.599        0.000
## 22    H68      0.000      0.001        0.397        0.521        0.000
## 23    H71      0.000      0.036        0.796        0.885        0.000
## 24     H9      0.000      0.002        0.348        0.466        0.000
##    carinulata_H carinata_L carinata_H elongata sublineata carinulata carinata
## 1         0.001      0.489      0.620    0.000      0.445      0.000    0.555
## 2         0.001      0.130      0.247    0.028      0.785      0.000    0.187
## 3         0.001      0.071      0.156    0.000      0.888      0.000    0.111
## 4         0.001      0.260      0.381    0.001      0.680      0.000    0.320
## 5         0.001      0.088      0.182    0.000      0.867      0.000    0.133
## 6         0.001      0.296      0.428    0.000      0.639      0.000    0.360
## 7         0.001      0.762      0.856    0.034      0.155      0.000    0.811
## 8         0.001      0.607      0.729    0.000      0.331      0.000    0.668
## 9         0.001      0.223      0.330    0.004      0.720      0.000    0.276
## 10        0.001      0.103      0.185    0.000      0.858      0.000    0.142
## 11        0.004      0.125      0.235    0.003      0.817      0.000    0.179
## 12        0.001      0.737      0.846    0.000      0.207      0.000    0.793
## 13        0.001      0.198      0.357    0.000      0.725      0.000    0.274
## 14        0.001      0.180      0.279    0.000      0.772      0.000    0.228
## 15        0.056      0.655      0.777    0.050      0.198      0.035    0.717
## 16        0.001      0.172      0.272    0.000      0.779      0.000    0.221
## 17        0.001      0.113      0.202    0.000      0.844      0.000    0.156
## 18        0.001      0.496      0.621    0.001      0.440      0.000    0.559
## 19        0.001      0.100      0.189    0.020      0.837      0.000    0.143
## 20        0.001      0.588      0.709    0.000      0.350      0.000    0.650
## 21        0.001      0.401      0.527    0.000      0.536      0.000    0.464
## 22        0.001      0.479      0.603    0.000      0.459      0.000    0.541
## 23        0.001      0.096      0.190    0.017      0.843      0.000    0.140
## 24        0.001      0.533      0.651    0.000      0.407      0.000    0.593
##    popID
## 1   28NM
## 2   26NM
## 3   26NM
## 4   26NM
## 5   26NM
## 6   27NM
## 7   15TX
## 8   20NM
## 9   20NM
## 10  28NM
## 11  20NM
## 12  16TX
## 13  21NM
## 14  26NM
## 15  28NM
## 16  26NM
## 17  26NM
## 18  26NM
## 19  26NM
## 20  18TX
## 21  28NM
## 22  28NM
## 23  26NM
## 24  28NM
car_sub_hyb %>%
  summarise(mean(sublineata), nrow(.))
##   mean(sublineata) nrow(.)
## 1        0.6075833      24
car_elong_hyb <- ci %>%
  dplyr::filter(carinata_L > t_low & elongata_L> t_low & sublineata_L < t_low) %>%
  arrange(desc(elongata_L))
car_elong_hyb
##    sampID elongata_L elongata_H sublineata_L sublineata_H carinulata_L
## 1     E29      0.221      0.305            0        0.002            0
## 2     G87      0.082      0.137            0        0.005            0
## 3 G01_rep      0.078      0.129            0        0.001            0
##   carinulata_H carinata_L carinata_H elongata sublineata carinulata carinata
## 1        0.001      0.694      0.779    0.262      0.000          0    0.737
## 2        0.001      0.862      0.917    0.109      0.001          0    0.890
## 3        0.001      0.870      0.922    0.103      0.000          0    0.897
##   popID
## 1  15TX
## 2  15TX
## 3  15TX
car_elong_hyb %>%
  summarise(mean(carinata), nrow(.))
##   mean(carinata) nrow(.)
## 1      0.8413333       3
elong_sub_hyb <- ci %>%
  dplyr::filter(carinata_L < t_low  & sublineata_L > t_low & elongata_L > t_low)
elong_sub_hyb
##   sampID elongata_L elongata_H sublineata_L sublineata_H carinulata_L
## 1    G10      0.088      0.133        0.866        0.912            0
## 2    G50      0.092      0.140        0.860        0.908            0
## 3    H26      0.198      0.263        0.736        0.802            0
##   carinulata_H carinata_L carinata_H elongata sublineata carinulata carinata
## 1        0.001          0      0.003    0.110      0.890          0    0.001
## 2        0.001          0      0.001    0.115      0.885          0    0.000
## 3        0.001          0      0.001    0.230      0.770          0    0.000
##   popID
## 1  26NM
## 2  20NM
## 3  21NM
elong_sub_hyb %>%
  summarise(mean(sublineata), nrow(.))
##   mean(sublineata) nrow(.)
## 1        0.8483333       3

Vizualize bi-parental crosses!

This could be better represented as a simplex (triangle)?

pCS <- ggplot(car_sub_hyb) +
  geom_histogram(aes(carinata),
                 alpha = .5, fill = legend[which(legend$spp=="carinata"), "cbbPalette"], stat = "bin", binwidth = .05) +
  geom_histogram(aes(sublineata),
                 alpha = .5, fill = legend[which(legend$spp=="sublineata"), "cbbPalette"], stat = "bin", binwidth = .05) +
  ylim(c(0,8)) +
  scale_x_continuous(name= expression( ~italic("D. carinata")~ "X"  ~italic("D. sublineata")~ "q-values")) +
  theme_bw()
pCS

pCE <-ggplot(car_elong_hyb) +
  geom_histogram(aes(carinata),
                 alpha = .5, fill = legend[which(legend$spp=="carinata"), "cbbPalette"], stat = "bin", binwidth = .05) +
  geom_histogram(aes(elongata),
                 alpha = .5, fill = legend[which(legend$spp=="elongata"), "cbbPalette"], stat = "bin", binwidth = .05) +
  ylim(c(0,8)) +
    scale_x_continuous(name= expression( ~italic("D. elongata")~ "X"  ~italic(" D. carinata")~ "q-values")) +
  theme_bw()
pCE

pES <-ggplot(elong_sub_hyb) +
  geom_histogram(aes(elongata),
                 alpha = .5,fill = legend[which(legend$spp=="elongata"), "cbbPalette"], stat = "bin", binwidth = .05) +
  geom_histogram(aes(sublineata),
                 alpha = .5, fill = legend[which(legend$spp=="sublineata"), "cbbPalette"], stat = "bin", binwidth = .05) +
  ylim(c(0,8)) +
      scale_x_continuous(name= expression( ~italic("D. elongata")~ "X"  ~italic(" D. sublineata")~ "q-values")) +
  theme_bw()

cowplot::plot_grid(pCS, pCE, pES, nrow = 1, labels = c('A','B','C'))

ggsave("hybrid_distributions.jpg", dpi = 300, height = 4, width = 10, units = "in")

Which sites have individuals with admxiture?

sample_size <- samplesize %>%
  select('popID', 'n')

hybrids <- rbind(elong_sub_hyb, car_elong_hyb, carinu_hyb, car_sub_hyb)
nrow(hybrids)
## [1] 32
no_hybrids <- anti_join(ci, hybrids)
## Joining, by = c("sampID", "elongata_L", "elongata_H", "sublineata_L", "sublineata_H", "carinulata_L", "carinulata_H", "carinata_L", "carinata_H", "elongata", "sublineata", "carinulata", "carinata", "popID")
pure_carinu <- dplyr::filter(no_hybrids, carinulata > .9)  
 nrow(pure_carinu)
## [1] 179
pure_sub <- dplyr::filter(no_hybrids, sublineata > .9)  
   nrow(pure_sub)
## [1] 176
pure_carinata <- dplyr::filter(no_hybrids, carinata > .9 )  
   nrow(pure_carinata)
## [1] 69
pure_elongata <- dplyr::filter(no_hybrids, elongata > .9 )  
   nrow(pure_elongata)   
## [1] 93
pure_parentals <- rbind(pure_carinu, pure_sub, pure_carinata, pure_elongata)

what_else <- anti_join( ci, pure_parentals) %>%
  anti_join(., hybrids)
## Joining, by = c("sampID", "elongata_L", "elongata_H", "sublineata_L", "sublineata_H", "carinulata_L", "carinulata_H", "carinata_L", "carinata_H", "elongata", "sublineata", "carinulata", "carinata", "popID")
## Joining, by = c("sampID", "elongata_L", "elongata_H", "sublineata_L", "sublineata_H", "carinulata_L", "carinulata_H", "carinata_L", "carinata_H", "elongata", "sublineata", "carinulata", "carinata", "popID")
 num_admixed_inds <- hybrids %>%
  group_by(popID) %>%
  count(popID)
colnames(num_admixed_inds)[2] <- 'num_admixed'

num_admixed_inds
## # A tibble: 9 x 2
## # Groups:   popID [9]
##   popID num_admixed
##   <chr>       <int>
## 1 15TX            4
## 2 16TX            1
## 3 18TX            1
## 4 20NM            4
## 5 21NM            2
## 6 26NM           11
## 7 27NM            1
## 8 28NM            7
## 9 6KS             1
merge(sample_size, num_admixed_inds) %>%
  mutate(perc_admixed = num_admixed/n) %>%
  arrange(desc(perc_admixed))
##   popID  n num_admixed perc_admixed
## 1  26NM 18          11   0.61111111
## 2  28NM 17           7   0.41176471
## 3  15TX 16           4   0.25000000
## 4  20NM 17           4   0.23529412
## 5  21NM 14           2   0.14285714
## 6   6KS 10           1   0.10000000
## 7  18TX 17           1   0.05882353
## 8  16TX 18           1   0.05555556
## 9  27NM 20           1   0.05000000
hybrids %>%
  filter_at(vars(elongata:carinata),
            ~ . < .95) %>%
  group_by(popID) %>%
  summarize_at(vars(carinata:sublineata), funs(mean, var))
## Warning: `funs()` was deprecated in dplyr 0.8.0.
## Please use a list of either functions or lambdas: 
## 
##   # Simple named list: 
##   list(mean = mean, median = median)
## 
##   # Auto named with `tibble::lst()`: 
##   tibble::lst(mean, median)
## 
##   # Using lambdas
##   list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
## # A tibble: 9 x 7
##   popID carinata_mean carinulata_mean sublineata_mean carinata_var
## * <chr>         <dbl>           <dbl>           <dbl>        <dbl>
## 1 15TX          0.834          0                0.039      0.00568
## 2 16TX          0.793          0                0.207     NA      
## 3 18TX          0.65           0                0.35      NA      
## 4 20NM          0.281          0                0.688      0.0797 
## 5 21NM          0.137          0                0.748      0.0375 
## 6 26NM          0.200          0                0.784      0.0205 
## 7 27NM          0.36           0                0.639     NA      
## 8 28NM          0.468          0.0999           0.425      0.0397 
## 9 6KS           0              0.088            0.912     NA      
## # … with 2 more variables: carinulata_var <dbl>, sublineata_var <dbl>

Which sites have multiple spp, even without admixture?

foo <- q_popmap %>%
  filter_all(any_vars(. > .95)) %>%
  group_by(popID) %>%
  summarize_at(vars(elongata:carinata), funs(sum)) %>%
  merge(sample_size)

max_ancestry <-  apply(X=foo[,2:5], MARGIN=1, FUN=max)
### if a cluster has less than foo$sample_size-.5, will return any pop with an ind less than 50% of majority ancestry
mixed_pops <- foo[foo$n - .5 > max_ancestry,]
mixed_pops
##    popID elongata sublineata carinulata carinata  n
## 3   15TX    1.580      0.158      0.005   14.252 16
## 4   16TX   16.996      0.209      0.000    0.794 18
## 5   18TX    0.002     14.249      2.092    0.655 17
## 8   20NM    0.125     15.647      0.003    1.223 17
## 9   21NM    0.328     13.295      0.019    0.355 14
## 10  26NM    0.338     15.372      0.019    2.270 18
## 11  27NM    0.321     19.176      0.002    0.495 20
## 12  28NM    0.094      3.035      8.596    5.274 17
## 32   6KS    0.000      0.915      1.088    7.995 10
perc_admixed <- merge(sample_size, num_admixed_inds) %>%
  mutate(perc_admixed = num_admixed/n) %>%
  arrange(desc(perc_admixed))
perc_admixed
##   popID  n num_admixed perc_admixed
## 1  26NM 18          11   0.61111111
## 2  28NM 17           7   0.41176471
## 3  15TX 16           4   0.25000000
## 4  20NM 17           4   0.23529412
## 5  21NM 14           2   0.14285714
## 6   6KS 10           1   0.10000000
## 7  18TX 17           1   0.05882353
## 8  16TX 18           1   0.05555556
## 9  27NM 20           1   0.05000000
mix_admix <- full_join(mixed_pops, perc_admixed[,c('popID', 'perc_admixed')]) %>%
  arrange(desc(perc_admixed))
## Joining, by = "popID"
mix_admix
##   popID elongata sublineata carinulata carinata  n perc_admixed
## 1  26NM    0.338     15.372      0.019    2.270 18   0.61111111
## 2  28NM    0.094      3.035      8.596    5.274 17   0.41176471
## 3  15TX    1.580      0.158      0.005   14.252 16   0.25000000
## 4  20NM    0.125     15.647      0.003    1.223 17   0.23529412
## 5  21NM    0.328     13.295      0.019    0.355 14   0.14285714
## 6   6KS    0.000      0.915      1.088    7.995 10   0.10000000
## 7  18TX    0.002     14.249      2.092    0.655 17   0.05882353
## 8  16TX   16.996      0.209      0.000    0.794 18   0.05555556
## 9  27NM    0.321     19.176      0.002    0.495 20   0.05000000

How many 'pure' sites did we detect?

setdiff(popmap$popID, mix_admix$popID )
##  [1] "11CO"       "12AZ"       "19TX"       "1WY"        "2CO"       
##  [6] "31NV"       "32UT"       "33NV"       "34UT"       "37CR"      
## [11] "38CR"       "39CR"       "41CR"       "43GR"       "44GR"      
## [16] "46CH"       "47TX"       "49TX"       "4CO"        "50TX"      
## [21] "51TX"       "52TX"       "5CO"        "8OK"        "CARINA_LAB"
## [26] "SUB_LAB"
length(setdiff(popmap$popID, mix_admix$popID ))
## [1] 26

Pull pure parental species for reference down the road.

pure_parentals <- ci %>%
  dplyr::filter(carinata_L > 1 - t_low | elongata_L > 1 - t_low |
           sublineata_L > 1 - t_low | carinulata_L > 1 - t_low) %>%
  merge(popmap)

pure_parentals$assignment <- colnames(pure_parentals[,c('elongata', 'carinata', 'carinulata', 'sublineata')])[max.col(pure_parentals[,c('elongata', 'carinata', 'carinulata', 'sublineata')])]

table(pure_parentals$assignment)
## 
##   carinata carinulata   elongata sublineata 
##         64        179         92        166
pure_parentals_popmap <- pure_parentals[,c('sampID', 'assignment')]
write.table(pure_parentals_popmap, "../info/pure_parentals_popmap.tsv", col.names = F, quote = F, row.names = F, sep = "\t")

Substructure

Subset individuals for substructure of elongata, carinulata, and everything but carinulata ("others").

elongata <- dplyr::filter(ci, elongata_L > t_low) %>%
  select('sampID', 'popID')
nrow(elongata)
## [1] 99
write.table(elongata,
  file = "../info/elongata_popmap.tsv",
            sep = "\t",
            quote=F, col.names = F, row.names = F)

carinulata <- dplyr::filter(ci, carinulata_L > t_low) %>%
  select('sampID', 'popID')
nrow(carinulata)
## [1] 182
write.table(carinulata,
  file = "../info/carinulata_popmap.tsv",
            sep = "\t",
            quote=F, col.names = F, row.names = F)

others_popmap <- anti_join(popmap, carinulata)
## Joining, by = c("sampID", "popID")
write.table(others_popmap,
  file = "../info/others_popmap.tsv",
            sep = "\t",
            quote=F, col.names = F, row.names = F)

Population structure within D. carinulata

sfiles <- list.files(path="../substructure/carinulata/structure_analysis/Results/",
  pattern=glob2rx("outfile_*f"),full.names=TRUE)

carinu_slist <- readQ(sfiles,filetype="structure",
               indlabfromfile = TRUE,
               readci = TRUE)

popmap <- left_join(data.frame(sampID = rownames(carinu_slist[[1]])),
                    masterpopmap)
## Joining, by = "sampID"
popID <- data.frame(popID = popmap$popID)
popID <- data.frame(lapply(popID, as.character), stringsAsFactors=FALSE)
sapply(popID, is.character)
## popID 
##  TRUE
pop_order <- sample_meta.dat[,c("popID", "Order")]
pop_order <- semi_join(pop_order, popID) %>%
  arrange(Order)
## Joining, by = "popID"
pop_order <- pop_order[,1]

legend_colors
##   cbbPalette human_colors            spp
## 1    #56B4E9         blue       carinata
## 2    #E69F00       orange       elongata
## 3    #009E73        green     sublineata
## 4    #F0E442       yellow     carinulata
## 5    #0072B2    dark blue   elongata_sub
## 6    #D55E00          red carinulata_sub
## 7    #9370DB       purple          other
carinu_legend <- legend_colors[c(4,6,7),]
carinu_cbbPalette <- carinu_legend[,'cbbPalette']
carinu_cbbPalette
## [1] "#F0E442" "#D55E00" "#9370DB"
sr1 <- summariseQ(tabulateQ(carinu_slist))
p <- evannoMethodStructure(data=sr1,exportplot=F,returnplot=T,returndata=F,basesize=12,linesize=0.7)
grid.arrange(p)

evannoMethodStructure(data=sr1,
                  exportplot=T,
                  exportpath = "../substructure/carinulata/structure_analysis/plots/",
                outputfilename = "carinulata_deltaK",
                returndata=F,
                height = 10, width = 9)
## ../substructure/carinulata/structure_analysis/plots//carinulata_deltaK.png exported.
plotQ(alignK(carinu_slist[c(21:30)]),
            imgoutput="join",
                  exportpath = "../substructure/carinulata/structure_analysis/plots/",
  outputfilename =  "Joined-K2",  
  exportplot=T,
  showyaxis=T,
  showsp = F,
  clustercol=carinu_cbbPalette,
  subsetgrp = pop_order,
  height = .1,
  ordergrp = T,
  grplab = popID,
  linepos = 1,
  divtype = 1,
  grplabangle=45,
  grplabpos = .8,
  grplabspacer = 0,
  grplabheight=70,
  width=100,
  showtitle=T,
  grplabsize = 4,
  titlelab="D. carinulata Population Structure, K=2",
  titlesize = 20)
## Drawing plot ...
## ../substructure/carinulata/structure_analysis/plots//Joined-K2.png exported.
plotQ(alignK(carinu_slist[c(31:40)]),
            imgoutput="join",
                  exportpath = "../substructure/carinulata/structure_analysis/plots/",
  outputfilename =  "Joined-K3",  
  exportplot=T,
  showyaxis=T,
  showsp = F,
  clustercol=carinu_cbbPalette,
  subsetgrp = pop_order,
  height = .1,
  ordergrp = T,
  grplab = popID,
  linepos = 1,
  divtype = 1,
  grplabangle=45,
  grplabpos = .8,
  grplabspacer = 0,
  grplabheight=70,
  width=100,
  showtitle=T,
  grplabsize = 4,
  titlelab="D. carinulata Population Structure, K=3",
  titlesize = 20)
## Drawing plot ...
## ../substructure/carinulata/structure_analysis/plots//Joined-K3.png exported.
k4 <- plotQ(alignK(carinu_slist[c(41:50)]),
            imgoutput="join",
            returnplot=T,exportplot=F,basesize=11,
            clustercol=cbbPalette,
            subsetgrp = pop_order,
            ordergrp = T,
            grplab = popID,
            grplabangle=-90)
grid.arrange(k4$plot[[1]])

k5 <- plotQ(alignK(carinu_slist[c(51:60)]),
            imgoutput="join",
            subsetgrp = pop_order,
            returnplot=T,exportplot=F,basesize=11,
            clustercol=cbbPalette,
            ordergrp = T,
            grplab = popID,
            grplabangle=-90)
grid.arrange(k5$plot[[1]])

k6 <- plotQ(alignK(carinu_slist[c(61:70)]),
            imgoutput="join",
            subsetgrp = pop_order,
            returnplot=T,exportplot=F,basesize=11,
            clustercol=cbbPalette,
            ordergrp = T,
            grplab = popID,
            grplabangle=-90)
grid.arrange(k6$plot[[1]])

k7 <- plotQ(alignK(carinu_slist[c(71:80)]),
            imgoutput="join",
            subsetgrp = pop_order,
            returnplot=T,exportplot=F,basesize=11,
            clustercol=cbbPalette,
            ordergrp = T,
            grplab = popID,
            grplabangle=-90)
grid.arrange(k7$plot[[1]])

plotQ(carinu_slist[32],
  exportpath = "../substructure/carinulata/structure_analysis/plots/",
  outputfilename =  "MainTextFigure_carinuK3",  
  exportplot=T,
  showyaxis=T,
  clustercol=carinu_cbbPalette,
  subsetgrp = pop_order,
  sortind="all",
  grplab = popID,
  linepos = 1,
  linesize=0.8,
  pointsize = 5,
  divtype = 1,
  divsize = 1.2,
  grplabangle= 90,
  grplabpos = .9,
  grplabspacer = 0,
  grplabheight=40,
  grplabsize = 25,
  height = .05,
  width=100,
  showlegend = T,
  showsp = F, 
  legendkeysize=40,legendtextsize=40, legendpos="left", legendlab= c("D. carinulata Fukang Ecotype", "D. carinulata Chilik Ecotype", "Other"),
  basesize=40)
## Drawing plot ...
## ../substructure/carinulata/structure_analysis/plots//MainTextFigure_carinuK3.png exported.

Population structure within D. elongata

sfiles <- list.files(path="../substructure/elongata/structure_analysis/Results/",
  pattern=glob2rx("outfile_*f"),full.names=TRUE)

elong_slist <- readQ(sfiles,filetype="structure",
               indlabfromfile = TRUE,
               readci = TRUE)

popmap <- left_join(data.frame(sampID = rownames(elong_slist[[1]])),
                    masterpopmap)
## Joining, by = "sampID"
popID <- data.frame(popID = popmap$popID)
popID <- data.frame(lapply(popID, as.character), stringsAsFactors=FALSE)
sapply(popID, is.character)
## popID 
##  TRUE
pop_order <- sample_meta.dat[,c("popID", "Order")]
pop_order <- semi_join(pop_order, popID) %>%
  arrange(Order)
## Joining, by = "popID"
pop_order <- pop_order[,1]

legend_colors
##   cbbPalette human_colors            spp
## 1    #56B4E9         blue       carinata
## 2    #E69F00       orange       elongata
## 3    #009E73        green     sublineata
## 4    #F0E442       yellow     carinulata
## 5    #0072B2    dark blue   elongata_sub
## 6    #D55E00          red carinulata_sub
## 7    #9370DB       purple          other
elong_legend <- legend_colors[c(2,5,7),]
elong_cbbPalette <- elong_legend[,'cbbPalette']
elong_cbbPalette
## [1] "#E69F00" "#0072B2" "#9370DB"
sr1 <- summariseQ(tabulateQ(elong_slist))

evannoMethodStructure(data=sr1,
                  exportplot=T,
                  exportpath = "../substructure/elongata/structure_analysis/plots/",
                outputfilename = "elongata_deltaK",
                returndata=F,
                height = 10, width = 9)
## ../substructure/elongata/structure_analysis/plots//elongata_deltaK.png exported.
k2 <- plotQ(alignK(elong_slist[c(21:30)]),
            imgoutput="join",
            returnplot=T,exportplot=F,basesize=11,
            clustercol=cbbPalette,
            subsetgrp = pop_order,
            ordergrp = T,
            grplab = popID,
            linepos = 1,
            grplabangle=-45,
            grplabpos = 0,
            grplabheight=20)
grid.arrange(k2$plot[[1]])

k3 <- plotQ(alignK(elong_slist[c(31:40)]),
            imgoutput="join",
            returnplot=T,exportplot=F,basesize=11,
            clustercol=cbbPalette,
            ordergrp = T,
            grplab = popID,
            grplabangle=-90)
grid.arrange(k3$plot[[1]])

k4 <- plotQ(alignK(elong_slist[c(41:50)]),
            imgoutput="join",
            returnplot=T,exportplot=F,basesize=11,
            clustercol=cbbPalette,
            ordergrp = T,
            grplab = popID,
            grplabangle=-90)
grid.arrange(k4$plot[[1]])

k5 <- plotQ(alignK(elong_slist[c(51:60)]),
            imgoutput="join",
            returnplot=T,exportplot=F,basesize=11,
            clustercol=cbbPalette,
            ordergrp = T,
            grplab = popID,
            grplabangle=-90)
grid.arrange(k5$plot[[1]])

k6 <- plotQ(alignK(elong_slist[c(61:70)]),
            imgoutput="join",
            returnplot=T,exportplot=F,basesize=11,
            clustercol=cbbPalette,
            ordergrp = T,
            grplab = popID,
            grplabangle=-90)
grid.arrange(k6$plot[[1]])

k7 <- plotQ(alignK(elong_slist[c(71:80)]),
            imgoutput="join",
            returnplot=T,exportplot=F,basesize=11,
            ordergrp = T,
            grplab = popID,
            grplabangle=-90)
grid.arrange(k7$plot[[1]])

plotQ(elong_slist[37],
            exportpath = "../substructure/elongata/structure_analysis/plots/",
        outputfilename =  "MainTextFigure_elongataK3",  
  exportplot=T,
  showyaxis=F,
  clustercol=elong_cbbPalette,
  subsetgrp = pop_order,
  sortind="all",
  grplab = popID,
  linepos = 1,
  linesize=0.8,
  pointsize = 5,
  divtype = 1,
  divsize = 1.2,
  grplabangle= 90,
  grplabpos = .9,
  grplabspacer = 0,
  grplabheight=40,
  grplabsize = 25,
  height = .05,
  width=100,
  showlegend = F,
  showsp = F,
  basesize=40)
## Drawing plot ...
## ../substructure/elongata/structure_analysis/plots//MainTextFigure_elongataK3.png exported.

Is there a latitudinal gradient found for ecotypes or ancestry?

library(ggtext) 

global_q.dat <- slist[[41]]
global_q.dat$sampID <- rownames(global_q.dat)
global_q.dat.popmap.meta <- merge(global_q.dat,
                    masterpopmap) %>%
  merge(sample_meta.dat)   %>%
 dplyr::filter(Grp == "potential_hybrid_zone")

sub_plot <- ggplot(global_q.dat.popmap.meta, aes(x = Latitude, y = Cluster2)) +
  geom_smooth(method = "lm", color = "red") +
  geom_jitter(aes(color = popID, shape = cat), alpha = .8, size = 4) +
  ylim(c(0,1)) +
  ylab("Proportion of *D. sublineata* ancestry") +
  guides(col = guide_legend(ncol = 2)) +
  theme_bw(base_size = 12)+
  theme(axis.title.y = element_markdown())
sub_plot
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 6 rows containing missing values (geom_smooth).
## Warning: Removed 83 rows containing missing values (geom_point).

ggsave("sublineata_ancestry_latitude.jpg")
## Saving 7 x 5 in image
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 6 rows containing missing values (geom_smooth).
## Warning: Removed 89 rows containing missing values (geom_point).
fit_carina <- lm(global_q.dat.popmap.meta$Cluster4 ~ global_q.dat.popmap.meta$Latitude)
summary(fit_carina)
## 
## Call:
## lm(formula = global_q.dat.popmap.meta$Cluster4 ~ global_q.dat.popmap.meta$Latitude)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.81015 -0.26457  0.06703  0.08528  0.67864 
## 
## Coefficients:
##                                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                       -3.042344   0.251300  -12.11   <2e-16 ***
## global_q.dat.popmap.meta$Latitude  0.101408   0.007729   13.12   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3187 on 275 degrees of freedom
## Multiple R-squared:  0.385,  Adjusted R-squared:  0.3827 
## F-statistic: 172.1 on 1 and 275 DF,  p-value: < 2.2e-16
fit_sub <- lm(global_q.dat.popmap.meta$Cluster2 ~ global_q.dat.popmap.meta$Latitude)
summary(fit_sub)
## 
## Call:
## lm(formula = global_q.dat.popmap.meta$Cluster2 ~ global_q.dat.popmap.meta$Latitude)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.75212 -0.12038 -0.01282  0.24788  0.94804 
## 
## Coefficients:
##                                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                        4.570464   0.270919   16.87   <2e-16 ***
## global_q.dat.popmap.meta$Latitude -0.121256   0.008333  -14.55   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3436 on 275 degrees of freedom
## Multiple R-squared:  0.435,  Adjusted R-squared:  0.433 
## F-statistic: 211.8 on 1 and 275 DF,  p-value: < 2.2e-16
plot(fit_sub, which=1, col=c("blue"))

plot(fit_sub, which=2, col=c("blue"))

global_q.dat.popmap.meta$predicted <- predict(fit_sub)   # Save the predicted values
global_q.dat.popmap.meta$residuals <- residuals(fit_sub) # Save the residual values

ggplot(global_q.dat.popmap.meta, aes(x = Latitude, y = Cluster2)) +
  geom_smooth(method = "lm", se = FALSE, color = "lightgrey") +     # regression line  
  geom_segment(aes(xend = Latitude, yend = predicted), alpha = .2) +      # draw line from point to line
  geom_point(aes(color = abs(residuals), size = abs(residuals))) +  # size of the points
  scale_color_continuous(low = "green", high = "red") +             # colour of the points mapped to residual size - green smaller, red larger
  guides(color = FALSE, size = FALSE) +                             # Size legend removed
  geom_point(aes(y = predicted), shape = 1) +
  ylab("Proportion of D. sublineata ancestry") +
  theme_bw()
## `geom_smooth()` using formula 'y ~ x'

ggsave("sublineata_ancestry_latitude_fit.jpg")
## Saving 7 x 5 in image
## `geom_smooth()` using formula 'y ~ x'
carinu_q.dat <- carinu_slist[[32]]
carinu_q.dat$sampID <- rownames(carinu_q.dat)
carinu_q.popmap.meta <- merge(carinu_q.dat,
                    masterpopmap) %>%
  merge(sample_meta.dat) #%>%
  # dplyr::filter(Cluster3 < .5)  # %>%
  # dplyr::filter( popID != "46CH"  & popID != "1WY"& popID != "18TX" &  popID != "28NM")


carinu_plot <- ggplot(carinu_q.popmap.meta) +
  geom_smooth(aes(Latitude, Cluster2), method = "lm", color = "red") +
  geom_jitter(aes(x = Latitude, y = Cluster2, color = popID, shape = cat), alpha = .8, size = 4) +
  ylab("Proportion of Chilik ancestry") +
  ylim(c(0,1)) +
  guides(col = guide_legend(ncol = 2)) +
  theme_bw(base_size = 12) 
carinu_plot
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 2 rows containing missing values (geom_smooth).
## Warning: Removed 1 rows containing missing values (geom_point).

ggsave("fukang_ancestry_latitude.jpg")
## Saving 7 x 5 in image
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 2 rows containing missing values (geom_smooth).
fit_carinu <- lm(carinu_q.popmap.meta$Cluster2 ~ carinu_q.popmap.meta$Latitude)
summary(fit_carinu)
## 
## Call:
## lm(formula = carinu_q.popmap.meta$Cluster2 ~ carinu_q.popmap.meta$Latitude)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.86712 -0.47929  0.02838  0.40506  0.56583 
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    2.641788   0.297754   8.872 7.03e-16 ***
## carinu_q.popmap.meta$Latitude -0.056325   0.007543  -7.468 3.39e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4096 on 180 degrees of freedom
## Multiple R-squared:  0.2365, Adjusted R-squared:  0.2323 
## F-statistic: 55.77 on 1 and 180 DF,  p-value: 3.394e-12
plot(fit_carinu, which=1, col=c("blue"))

plot(fit_carinu, which=2, col=c("blue"))

carinu_q.popmap.meta$predicted <- predict(fit_carinu)   # Save the predicted values
carinu_q.popmap.meta$residuals <- residuals(fit_carinu) # Save the residual values

ggplot(carinu_q.popmap.meta, aes(x = Latitude, y = Cluster2)) +
  geom_smooth(method = "lm", se = FALSE, color = "lightgrey") +     # regression line  
  geom_segment(aes(xend = Latitude, yend = predicted), alpha = .2) +      # draw line from point to line
  geom_point(aes(color = abs(residuals), size = abs(residuals))) +  # size of the points
  scale_color_continuous(low = "green", high = "red") +
  ylab("Proportion of Fukang ancestry") + # colour of the points mapped to residual size - green smaller, red larger
  guides(color = FALSE, size = FALSE) +      
    ylab("Proportion of Chilik ancestry") +# Size legend removed
  geom_point(aes(y = predicted), shape = 1) +
  theme_bw()
## `geom_smooth()` using formula 'y ~ x'

ggsave("fukang_ancestry_latitude_fit.jpg")
## Saving 7 x 5 in image
## `geom_smooth()` using formula 'y ~ x'
library(cowplot)
plots <- align_plots( sub_plot, carinu_plot, align = 'v', axis = 'l')
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 6 rows containing missing values (geom_smooth).
## Warning: Removed 86 rows containing missing values (geom_point).
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 2 rows containing missing values (geom_smooth).
plot_grid(plots[[1]],
          plots[[2]],
          ncol = 1, 
          labels = c('A', 'B'))

ggsave("sublineata_carinulata_ancestry_latitude.png", dpi = 300, height = 8, width = 8, units = "in")